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Crater formation: can macroscopic scaling laws be used in microscopic cratering? 
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Using classical molecular dynamics simulations we examine the formation of craters during 0.4 - 
100 keV Xe bombardment of Au. Our simulation results, and comparison with experiments and 
simulations of other groups, are used to examine to what extent analytical models can be used to 
predict the size and properties of craters. We do not obtain a fully predictive analytical model 
(with no fitting parameters) for the cratering probability, because of the difficulty in predicting the 
probability of cascades splitting into subcascades, and the relation of the heat spike lifetime and 
energy density. We do, however, demonstrate that the dependence of the crater size on the incident 
ion energy can be well understood qualitatively in terms of the lifetime of the heat spike and the 
cohesive energy of the material. We also show that a simple energy density criterion can not be 
used to predict cratering in a wide ion energy range because of the important role of the heat spike 
lifetime in high-energy cascades. The cohesive energy dependence differs from that obtained for 
macroscopic cratering (observed e.g. in astrophysics) because of the crucial role of melting in the 
development of heat spikes. 
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I. INTRODUCTION 

Surface modification of materials by incident ions 
has been observed using electron microscopy, scanning 
tunneling microscopy (STM) and atomic forces mi- 
croscopy (AFM)Era. |-.A large variety of features 
been studipcL hillocksO, depressionsacl, £rater rimst 
adatomsErEJ, and surface rougheningtJ. Hillocks can 
appear when an energetic process occurs a few layers 
below the surface. For instance, an energetic recoil can 
create a mini-spike which melts the surrounding region, 
creating a low deasity region of larger volume which 
raises the surfacelij. When the energy loss per unit 
path length of the projectile, dE/dx, and the sputtering 
yield are relatively small, adatoms are observed in 
both experiments and simulations. For larger energy 
deposition (and larger yields) a crater is formed. For 
even larger yields re-deposition of the ejecta plus plastic 
deformation occurs producing craters with runs, studied 
recently for ion bombardment of polymerscl. Craters 
are also produced by cluster ion bombardment in. which 
non-linear effects lead to enhanced sputteringi3~E2l. This 
has been studied in the velocity regime in which nuclear 
(elastic) energy loss dominates over electronic energy loss 
and has also been observed in experimentsEH lEI. 

Although several studies have used molecular dy- 
namics (MJ^pcp mp i utBr . simulations to examine crater 
formatiorilJ'ES'QHaEilfia^Ea, most of them have been 
limited to fairly low total ion or cluster impact energies, 
< 10 keV. In this energy range, a fairly good description 
of the mechanisms of cratering has been found in dense 
fee metalsoO. Yet MD simulations of mixing in the bulk 
have shown that the amount of atom displacements keeps 
increasing superlinearly up to ion energies of ~ 100 keV 
in dense fee metals, and that the Aeat spikes formed 
can persist for tens of picosecondsEJ. This poses the 
question whether the long lifetime of high-energy spikes 
can alter the mechanism of cratering in fee metals, as 
significant liquid flow might take place in cascades with 



long lifetimes. In fact, in a recent paper Aderjan and 
UrbassekB suggest that liquid flow might explain why the 
cratering cohesive energy dependencies differ from those 
predicted by macroscopic modelsEiTEil, but they do not 
give a description of the mechanism by which this might 



In this paper, we use MD simulations to examine 
crater formation by Xe recoils impacting on Au in an 
unprecedentedly wide energy range, ranging from 0.4 up 
to 100 keV of initial ion energy. We chose this ion-solid 
combination because TEM experiments are available for 
cratering in the same systemEjOEj. We find that at 
higher energies the fundamental mechanism leading to 
cratering does indeed change. While in the low-energy 
regime (< 10 keV) the mechanism can be understood 
on the basis of a high kinetic energy density alone, 
in agreement with previous models for dense metals, 
we show that at high energies (> 50 keV) cratering 
can result from lower kinetic energy densities due to 
the long lifetime of the heat spike. We also present 
analytical models for both energy regimes, and compare 
our simulation results with experiments and results in 
other types of materials. 

This paper is organized as follows. First some details of 
the simulation and crater identification and measurement 
are given in section |fi|. In section III results for 



a few individual events are presented, leading to the 
identification of weak points in current analytical models 
of cratering. In the following section we discuss the 
formation mechanisms of our craters in detail, and relate 
them to previous models. In section |v| we compare 
our results directly with experiments, and in section VI 
discuss why macroscopic scaling laws do not apply to 
atomic systems. We then compare our results for single- 
ion bombardment with results for cluster bombardment 



in section VII and finally in section VIII show that with 
appropriate scaling crater sizes in both metals and some 
organic solids can be understood in the same framework. 



II. METHOD 
A. Molecular dynamics simulations 

The basic MD simulation methods used in xkiw, pyjtfk 
has been described in several previous papers! 1 ffil'Ej'Ea, 
so in here we only recall the basic principles, and the 
features which differ from those. In simulating ion 
irradiation of a surface, we place an incident ion on 
a random position a few A above the surface, and 
give it a kinetic energy of 0.4 - 100 keV towards the 
sample. The incident angle is chosen in an off-channeling 
direction close to the surface normal. Most simulations 
were carried out for a (001) surface, but six 100 keV 
events were also simulated for a (111) surface. No major 
difference were observed in the crater sizes for the two 
surfaces for 100 keV energies. The development of the 
system of atoms is followed until the cascade has cooled 
down close to the ambient temperature, which for the 
cascades presented here was always K to definitely rule 
out any post-cascade damage annealing. A few 50 keV 
events (the results of which are not presented here) were 
also simulated at 300 K, and found to give similar crater 
sizes as the K events. Varying the initial position 
of the recoil atom can cause the resulting cascades to 
behave very differently depending on where the strongest 
collisions occur. 

The simulation cells had periodic boundaries in the 
x and y dimensions, and a fixed bottom layer in the z 
direction. They were cooled down to K using Berendsen 
temperature control at the cell borders and a few layers 
at the bottom above the fixed layer. The simulation cells 
had at least 16 atoms per eV of incident ion energy, which 
was enough to prevent cell heating beyond the melting 
point or pressure wave reflection from the borders strong 
enough to affect the cascade outcome. To ensure that no 
artificial border effects occurred, runs were automatically 
stopped and restarted in a larger simulation cell if a 
recoiling atom with an energy higher than 20|-e,V (which 
is less than the threshold displacement energyEa) entered 
the temperature scaling region. 

Most simulations were carried out in Au modeled by 
the embedded atom method (EAM) potential of FoilesEl, 
smoothly joined to the universal repulsive interatomic 
potential of Ziegler, Biersack and LittmarkO at small 
interatomic separationslij. To probe the cohesive energy 
dependence of the crater size, we used artificial modifica- 
tions of this potential for different values of the cohesive 
energy. To be precise, in the EAM formalismlfJo we 
scaled the pair potential and the embedding energy of 
the potential by a factor / in the range 0.5 - 1.6, while 
leaving the electron density unmodified. To preserve the 
ballistic properties of the potential, however, the factor 
/ was scaled smoothly back to 1.0 at small interatomic 
separations r, in the same r range where the repulsive 
potential was onset. This scaling gives an interatomic 
potential with the same equilibrium lattice constant as 
the original potential, but with the cohesive energy Uq 
and elastic moduli scaled by the factor /. We further 
verified by simulations that the melting point T mc i t scaled 
highly accurately by 1//, as expected from Lindemann's 



law§. 

For reference purposes, we note that at K the binding 
energy of our EAM Au is U = 3.930 eV, the lattice 
constant a — 4.080, and hence the atomic density is n = 
0.0589 A -3 . The melting temperature of the simulated 
solid is 1110±20K. 



B. Finding non-channeling directions 



To enable the use of simple binary collision approx- 
imation programs such as TRIMcJ for quick estimates 
of cascade energy densities and penetration depths, it 
is important to use an incident angle in the simulation 
corresponding to a non-channeling direction. To find 
such a direction, we used the MDRANGE code, which 
is an ion range calculation code which accounts for the 
crystal structurec3. 

TABLE I. Simulated mean range (R) and straggle Sr 
values for 50 keV Xe bombardment of (001) Au surfaces at 
K, using no zero-point atom displacements. Using thermal 
atom displacements corresponding to 300 K gave almost as 
strong channeling effects, for instance giving a mean range 
R = 103 ±3 A for 9 = 10°, <j> = 20-30°. The ion range is here 
defined as the ion penetration depth. 6 is the angle between 
the initial ion velocity vector and the (001) surface normal 
("tilt" angle), and <j> the angle between the initial vector and 
the (100) direction in the surface plane ("twist" angle). All 
angles are given in degrees. A range of angles such as "0-360" 
means the angle was selected randomly in this range. 



(p R (A) Sr (A) 

1240 ± 50 780 

5 20-30 670 ± 20 480 

10 20-30 108 ±5 110 

15 20-30 81 ±3 90 

20 20-30 79 ±3 60 

25 20-30 69 ±3 60 

30 20-30 75 ±2 76 

25 0-10 164 ± 14 251 

25 10-20 73 ±2 70 

25 20-30 69 ±3 60 

25 30-40 74 ±3 60 




FIG. 1. ((lard preprint n.b.: figure will be in better resolution in final published paper)) Crater formation by a 100 keV Xe 
ion hitting a (001) Au surface. The figure shows a cross-sectional slice of 8 atom layers in the (210) plane in the central part 
of the simulation cell. This cross-section was chosen to give a good illustration of the subcascade splitting. The arrow in the 
first frame indicates the initial direction of the incoming ion. The tip of the arrow shows the impact point of the ion projected 
on this cross section; the actual impact did not occur in these atom layers. The 0.25 ps and 1.0 ps snapshots show that the 
cascade splits into two almost separated subcascades during the ballistic phase of the cascade. The subcascade below the 
surface subsequently behaves like a cascade in tha-bulkEJ. The subcascade at the surface produces a crater between 2 and 40 
ps, and also causes a large atom cluster to sputteru. This sputtered cluster is so hot that it emits a large number of Au atoms 
and dimers, many of which redeposit on the surface. Notice also how an interstitial-like dislocation loop has formed at 20 rpe 
close to the right edge of the crater, and at 40 ps has produced an adatom island next to the crater by coherent displacementtJ. 
The final 150 ps snapshot shows the final crater structure, a vacancy loop produced by the subcascade inside the sample, and 
a complex dislocation structure on the right-hand side of the crater. 




FIG. 2. ((lanl preprint n.b.: figure will be in better 
resolution in final published paper)) (color) Crater produced 
in the event illustrated in Fig. hi seen from above. The colors 
indicate the height of the atoms. Blue and cyan atoms are 
below the original surface, and green, yellow and red atoms 
above, with the red atoms being highest up. The color scale 
has been chosen to emphasize atom layers at the surface; 
hence all atoms deeper than 1 unit cell (4.08 A) below the 
original surface have the same blue color. The maximum 
depth of the crater is about 44 A. Notice the regular adatom 
island (green atoms) on the lower right side of the crater, 
produced by the coherent displacement mechanism. The 
dislocation below the surface also produces a regular atom 
edge close to the crater, just next to the adatom island. The 
single adatoms far from the crater are atoms which have been 
redeposited on the surface from the sputtered atom clusters. 

In the particular case of heavy ion irradiation of (001) 
fee metals, we have observed that channeling effects can 
be extremely strong, requiring careful selection of both 
the tilt (9) and twist (</>) angles to obtain a good non- 
channeling direction. As we shall see below, this may 
be of crucial importance for the comparison of cratering 
probabilities with experiments, whence we give some 
sample range results for 50 keV Xe bombardment of (001) 
Au at K in Table |. The results show that the mean 
range R depends sensitively on both 9 and 4>. Ranges at 
other energies had a quite similar 9 and <j> dependence. 
To obtain an optimal non-channeling direction we used 
9 = (j) = 25° for the (001) surface at all energies. 



C. BCA Calculations 

A number of SRIM 2000.390 calculations were per- 
formed in order to compare binary collision approxima- 
tion (BCA) results with MD results. If good agreement 
is found e.g. in the energy densities, the SRIM code 
could be used to obtain quick estimates of cratering 
probabilities. The runs used the same incident angle, 9 = 
25°, U surf =2.6 eV, and U bu i k =3.9 eV. We used E D = 25 



eV as the displacement energy value for SRIMEl The 
low energy self-sputtering of Au at normal incidence 
was reproduced well with these values of the binding- 
energies and displacement energy. The total energy 
deposited up to certain depth, Eki n ,tot can be calculated 
by integrating the energy deposited as a function of 
depth z, Fd(z). Energy densities were calculated using a 
cylindrical volume (which contained the cascade better 
than a hemispherical volume). SRIM follows recoils 
until they reach the displacement energy, and therefore 
one should really compare with the MD energy density 
discarding atoms with Eki n > Ed- However, the energy 
deposition in SRIM represents the scenario before 0.1 ps, 
when there are only few energetic recoils and the cascade 
has not evolved significantly. 



D. Analysis 



Post-run analysis of the quenched MD simulation cells 
was performed to identify the surface features. If the last 
layer of atoms was located at z = 0, with the target in 
z < 0, all atoms in the range a/4 < z < r spu t were 
counted as adatoms, and all atoms above r sput were 
counted as sputtered atoms. r spu t was set to 20 A for 
all cases except the 100 keV ones, where r sput — 40 A. 
Both values are much larger than the cut-off radius of 
the potential, which is only 5.55 A. 

The crater and the crater rims were measured along 
two directions and the mean sizes were evaluated. The 
deviation of the crater shape from a circle was estimated 
with a parameter a = R < / R > , where i?< and i?> are the 
smaller/larger crater radii, respectively. This parameter 
was typically 0.6-0.9, since the craters were not circular 
but diamond shaped, following the symmetry of the 
lattice, as seen in cluster bombardment simulation of 
Cut! The area of the craters was therefore calculated as 
a circle and also as a parallelogram, which gives smaller 
areas, but within 10 - 30% of the value for a circular 
area. 

The rims were identified as structures with more 
than one atomic layer above the surface, to differentiate 
them from single-layer adatom islands and coherent 
displacement^, which occurs often at the side of the rim 
for the 20 - 100 keV events. 

The volume of the crater is approximated as the 
number of atoms 'excavated' from the target, nV = 
N a dat + N sp ut = N to t , where n is the equilibrium atomic 
density of the solid. If the number of atoms on top of 
the surface is used as an estimate of the crater volume, 
there are two small corrections which have not been 
taken into account. There may be atoms far from the 
crater which are on top of the surface because of coherent 
displacements. Besides, there may be interstitials or 
vacancies in the crater walls which change the normal 
density of the material, but their number was observed 
to be small, and in any case their effect seems to roughly 
cancel each other out. The 'measured' radius of the crater 
is not affected by these corrections. 



TABLE II. Results for individual cascade events simulated by MD and SRIM averages for the same energies. For MD, 
Ekin.tot is the total energy of the liquid atoms in the top 40 A, averaged over two times around 0.1 ps (for instance at 0.07 ps and 
0.15 ps). Ekin,mean is the mean energy of all atoms inside a cylinder of radius R ra d and length Depth. For the MD simulations 
Rlat, Rrad and "Depth" indicate the size of the cascade, e is the energy density calculated using a cylinder containing the 
cascade, calculated at 0.1 ps in the MD runs, n is the atomic density and Uo the equilibrium potential energy of the material. 
N vac is the number of vacancies, and Y is the initial sputtering yield. Note that because many of the sputtered atoms leave the 
surface in hot clusters, and can subsequently evaporate from the cluster and redeposit on the surface, the initial MD sputtering 
yield is not exactly comparable to experimentally measured yields. Individual simulations are labeled by an extra number, for 
instance, case number 8 for 10 keV bombardment is labeled 10-8. 





5 keV 


10 keV 


100 keV 


Event 


5 


5-2 


5-1 


5-4 


5-9 


10 


10-8 


100 


100-4 


Model 


SRIM 


MD 


MD 


MD 


MD 


SRIM 


MD 


SRIM 


MD 


Crater? 




no 


small 


yes 


yes 




yes 




yes 


Ekin.tot (keV) 


4.42 


2.125 


2.81 


2.71 


3.16 


8.32 


6.27 


25.4 


62.35 




1.84 


0.88 


1.17 


1.13 


1.31 


1.66 


2.4 


0.31 


0.76 


Rcrat 






8±3 


16±4 


16±4 




19.6±2 




36±6 


Rlat / Rrad 


13/18 


20 


25 


17 


12.5 


19/26 


20 


77/105 


60 


Depth 


30 


40 


30 


40 


25 


40 


35 


250 


60 


e/nUo 


0.8 


0.39 


0.43 


0.68 


2.35 


0.5 


0.62 


0.07 


0.21 


N vac 


119.5 


75 


73 


299 


152 


233.3 


408 


1,064.5 


4,538 


Y 


15.9 


2 


6 


31 


15 


21.3 


44 


36.3 


472 



TABLE III. Summary of simulation results for different energies. To calculate outflow time using Eq. |2| and spike times in 
Eq. |, we used R c = Ri for 1-20 keV, R c = 0.75ft for 50 keV, R c = 0.50ft for 100 keV. The spike times for the bulk of the 
MD simulation were calculated as the time when the number of liquid atoms was reduced by a factor of 10 from its maximum 
value, which occurs at « 0.5 ps. The spike times for the surface were calculated as the times when the number of liquid atoms 
at the surface had a maximum. Cratering probabilities were evaluated as explained in the text. 



E (keV) 


1 


5 


10 


20 


50 


100 


events 


20 


10 


20 


20 


10 


6 


ft (A)(MD) 


8.5±1.5 


13±4 


19.6±2 


24±4 


31±2.5 


33±2.5 


ft (A) (SRIM) 


7 


13 


19 


26 


43 


65 


t spike (ps)(MD)-bulk 


3.8 ±0.2 


9.7 ±0.4 


13 ±3 


13 ±1 


33±3 


30±10 


t spike (ps)(MD-surf)-no crater 


0.9 ±0.3 


2.3 ±0.4 


5 ±2 


18 ±2 


7.5 ±3 


4±2 


t ap ike (ps)(MD-surf)-crater 


0.9 ±0.2 


2.6 ±0.4 


6 ±2 


23 ±1 


12 ±3 


16±5 


t spike (ps)(Eq. B) 


0.8 


2.8 


6.0 


11.3 


17.3 


17.6 


touts (ps)(Eq. 2|) 


2.6 


4.9 


7.1 


9.8 


21.0 


21.2 


P [e (E ) > Ee] (MD) 


0.80 


0.78 


0.73 


0.7 


1 


1 


Psplit 


0.0 


0.1 


0.2 


0.3 


0.3 


0.4 


P (E a ) (MD) 


0.79 


0.67 


0.60 


0.55 


0.83 


0.60 


P (E ) (Eq. g) 


0.80 


0.70 


0.58 


0.49 


0.70 


0.60 



III. INDIVIDUAL CASCADES AND CRATERING 

A detailed picture of the cratering formation scenario 
can be obtained by considering individual MD simula- 
tions. For brevity, we use the following notation to denote 
individual events. The B'th event, were B is a running 
index, carried out with an energy of "A" keV is denoted 
by "A-B" . Thus for instance the eighth 10 keV event is 
marked by "10-8" . Results for a few events are listed in 
Table [n and discussed below. 

Fig. j] shows cross-sectional snapshots of a 100 keV 
cascade, event 100-4, with a crater being formed between 
2 and 40 ps. Fig. || shows the final crater produced 
in this event. The liquid flow of atoms builds the 
crater and crater rim. Notice that the cascade splits 
below the surface, and one of the subcascades does not 
reach the surface, even though it produces some coherent 
displacement^ of atoms next to the crater. Below we 
analyze in detail a few cascades in order to clarify some 
results regarding cratering probability, which is discussed 
in the following section. 

In Table || we compare MD results for individual 
events to average results from the SRIM 2000 BCA code. 
The MD values for cascade size, mean energy per atom, 
and so on compare reasonably well with SRIM, except 
for 100 keV. Note that the energy density at 0.1 ps, 
when the craters just start to form, is of the order of 
0.3 - 2uUq. For the 10 keV case shown in the table, 
by 1 ps this energy density decreased to ~ 0.075nE/o, 
in good agreement with other estimates of mean kinel 
energy (E^in) in the molten region of the cascadeo 
The case 10-8 forms an almost hemispherical crater, as 
well as the case 5-9. The case 5-4 has a large subcascade 
splitting, but both cascades are very close to the surface. 
On the other hand, case 5-2 has two subcascades, with 
some coherent atomic displacement on top of the larger 
cascade producing only few adatoms and a platelet. The 
case 100-4 forms a large crater, but notice that the energy 
density is lower than that for the 5 keV event that did 
not produce any crater. Since the cascade is so large, it 
can stay hot longer, and flow occurs during tens of ps, 
even though the energy density is low. Notice that the 
crater radius is much smaller than the cascade size, but 
the rim for 100-4 has a length of 100 ± 10 A, close to the 
lateral range of 105 A. 

From the analysis of the above cases and other events 
not discussed in detail here, we can already qualitatively 
conclude that: 

1. At low bombarding energies, cratering occurs for 
energy densities close to 0.5nUo. 

2. Energy density alone is not a good criterion for 
crater formation, since for large spikes the long 
lifetime of the spike will allow ejection even for 
energy densities much lower than nUo- To obtain a 
cratering probability one would need to take into 
account the probability of the spike being long 
lived, related to both spike radius and available 
energy. 

3. Cascade splitting should somehow be included, be- 
cause it can dilute the energy density significantly 



and concentrate it away from the surface. As a 
result the cratering probability is reduced. 



IV. FORMATION MECHANISMS 



Averaged simulation results can be found in Table 



III. The data for R c includes only the results of 
those simulations which produced a crater. The error 
bars indicate standard deviation of those events. The 
cratering probability was calculated as the ratio of 
cratering events to total number of events, P (E Q ) = 
N cra t/Ntotai, and it was always larger than 40% for the 
energies studied. 



A. Role of spikes at different energies 

In the early work of Thompson and Johai0, where 
large deviations from linear cascade theory were found, it 
was proposed that thermal spikes were not necessary, and 
that a decrease in the surface binding would be enough to 
explain the data. In our simulations spikes are certainly 
playing a large role and, since the surface damage is 
considerable, the surface binding also decreases. 

It has been claimed that the energy-. dens ity in the 
cascade determines the crater formationEnl3l23, and that 
an energy density larger than nUo is needed ptQ-,produce 
a crater. On the other hand, several modelscSfl assume 
that the flow of liquid atoms towards the surface can 
occur when atoms have Eki n ,mean = 3/csT m = 0.29 eV, 
giving a much lower energy density than before, e = 
0.07nlV We now use our simulation results to elucidate 
the reasons behind this apparent discrepancy. 

In order to verify the first assumption, a simple 
estimate can be made as follows. For instance, for E a = 
50 keV bombardment, we can assume an hemispherical 
crater of radius R c , and that E a is shared by all 
atoms inside the crater, which gain an energy E c , 
(2/3) nRlnE c = E a . Using R c = 30 A, close to the 
one found in MD (see Table |l|), gives N crat =3336, 
E c = 15eV/atom ~ 3.8C/o/atom, if most of the energy is 
originally deposited inside the crater region. N a d a toms ~ 
3370 from MD compares well with the missing crater 
volume. The agreement is quite good, except for energies 
below 20 keV, when craters are shallow, or energies 
above 50 keV, where a significant fraction of the energy 
is deposited deep inside the sample. However, as can 
be seen from Table f|, the energy density at energies 
above 10 keV is much smaller than uIIq, although higher 
than 0.07n[/o- At lower energies the energy density 
increases, but still is smaller than the minimum energy 
that would produce craters according to Ref. |l|. The 
cascade in the simulation 10-8, in Table y, has a crater 
size R cr at — 19.6 ± 2 A, and Eun,tot = 6.27 keV. Inside 
a cylinder of radius 20 A, and height 35 A there are 
2595 atoms, and Eki n ,mean — 2.4 eV < Uq. This gives 
e = 0.62?iC/ (nU = 0.2315 eV/A 3 ). Discarding the 
atoms with kinetic energy below Uq gives e c ~ 0.59nUo. 

Thus we see that the apparent contradiction in energy 
densities mentioned earlier may at least in part arise 



from the fact that the energy density needed for crater 
production does in fact strongly depend on the spike 
lifetime. 

Merkle and Jager proposed that surface spikes were the 
cause of the cratering and sputtering yield enhancement 
they observed in Bi + and Bi ++ bombardments of Au. 
Crater formation in our simulations is related to the 
probabil ity of the cascade being close to the surface (see 
Section [V C), as in their work. However, since they 
could only see craters with a radius larger than 25 A, 
they find a high energy threshold for crater formation. 
In our simulations we observe that crater formation can 
still occur at relatively low bombarding energies, but the 
crater size is well below the detection limit of Merkle 
and Jager (see Fig. ||). They also assume that atoms 
originally resident at the crater site are all sputtered, 
while we now know that most of the atoms will be 
redeposited at the crater rim and will not contribute to 
the sputtering yield. 

Based on the above discussion and the results in 



Section III, we can identify two regimes in the crater 
formation process. First, at low energy deposition only 
a relatively small hot region is created. For single ions 
this will occur when the penetration depth has roughly 
the same size as the lateral range of the ion, i.e. up 
to 10 - 20 keV in Xe — > Au. If the energy density 
deposited in the cascade is larger than 0.25nt7o a crater 
will be formed, provided the cascade is connected to the 
surface and did not decay into multiple subcascades. The 
probability of crater formation is slightly lower than one 
only because of the latter reason. For energies lower than 
1 keV, the lateral size of the cascade is of the order of the 
lattice spacing, and it would be difficult to distinguish a 
crater from a small vacancy cluster. We did perform a 
number of 400 eV events, and found that 4 out of 17 
events simulated resembled craters. However, because 
of the difficulty of clearly defining what is a crater with 
the very small number of atoms involved, we did not 
include these events in the quantitative analysis. In the 
low-energy (1-20 keV) regime, the crater has a roughly 
hemispherical shape, the radius of the crater is close to 
the radius of the lateral range of the ion in the solid, 
and the dependence of the size with energy follows the 

1 /3 

E law. Prompt sputtering of hot atoms cools down 
the spike significantly. 

Second, at high energy deposition, the spike region 
is large, and the energy density is lower than 0.25n[/o- 
However, since the spike is large, the center of the spike 
cools down much slower that its sides, and there are 
atoms which can flow out creating a crater. The crater 
radius is only a fraction of the lateral size of the cascade, 
since the borders cool down rapidly. Writing 



R c ^Ri [1.13-2.6 10- 5 (E o /U)] , 



(1) 




10 100 1000 
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FIG. 3. Crater radius R c as a function of the incident 
energy [Spr-i including experimental data of Donnelly and 
BirtcheiMy. 



B. Liquid Flow 

There will be a crater if the time for atomic outflow, 
tout/) is close to or larger than the lifetime of the spike, 
t spike- toutf can be estimated as the time it takes the 
atoms originally inside the crater to flow out. Assuming 
for simplicity a cylindrical outflow at constant velocity 
v c = y/2E c /m, and up to a depth R c , 

V C o n AR C 

n—7rR c t out f = mrR c =>• t outf = (2) 

4 V c 

v c is a mean velocity of the atoms after the outflow 

AR C 

begins. Writing E c = alio, toutf = 7=. = 

0.84 {R c /a) /yfa. (Ref. @). The spike lifetime can 
be estimated assuming lattice heat conductivity with a 
constant heat diffusivity k (Ref. B91). Then, 



t spiked ~ Rc/ (4«) 



(3) 



gives a very good estimate of the crater radius as a 
function of the lateral size of the cascade. 

Next we discuss this liquid flow in greater detail. 



If the lifetime of the spike is estimated as the time that 
it takes to cool down from the initial energy density to 
the critical energy density, t sp ik e can also be obtained 
analytically (assuming a Gaussian temperature profile), 
given the initial cascade parameters, and gives times 
within 30% of the simple estimate in Eq. H In order 
to check the validity of Eqs. || and ||, we compared 
the times from those equations to the MD results, using 
K, = 15 A 2 /ps (obtained from MD), and Eq. to obtain 
the spike radius. Note that typically it is assumed that 
r %ike x E (Ref. fl6|) but, following transport theory, the 
lateral range is roughly linear with E a which would mean 

T spike K Eo ■ 

Spike times are generally calculated neglecting the 
cooling due to evaporation. However, this cooling is 
significant when only a small hot region is created, and 
the penetration depth has roughly the same size as the 



lateral range of the ion, i.e. up to 10 - 20 keV in Xe 
— * Au. On the other hand, for 100 keV Xe — * Au, the 
lateral size of the cascade is ~ 60 A, while the ion range is 
- 135 A ( see Table [n]) . The flow of atoms to the surface 
starts at about 0.1 ps, even though the "ballistic" part 
is not finished. The ion is colliding further down from 
the surface, where the spike has not developed yet, but 
liquid flow already started above this region. 

As seen in Table [II, Eq. || understimates the 
spike lifetime in the bulk, because R c is used as the 
characteristic length scale. However, the agreement 
between Eq. |^ and the MD results for the lifetime 
of surface spikes in cratering events is good, while the 
surface spike lifetime in non-cratering events is shorter, 
especially at higher energies. The outflow times from 
Eq. U are close to the spike times from Eq. || for energies 
above 5 keV, making possible the crater formation by 
melt flow. 



C. Cratering probability 

The cratering probability for a projectile with energy 
E Q , P(E Q ), can be estimated as follows. 



P{E )= f P(E, 



, Ri) dRi 



where the probability of creating a crater of radius Ri 
is written as P (E , Ri). In order to simplify the analysis, 
we assume that the probability of producing a crater with 
certain radius Ri is a delta distribution for Ri = R c « Ri , 
P{E ,Ri) = P(E ,R c )S{Ri- R c ). This is especially 
valid at low energies, where fluctuations in crater size are 
not as large as for higher energy, -as was also noted in the 
experiment of Merkle and JagerEl. This probability can 
itself be split in three contributions. P [e (E Q ) > e c ] is the 
probability of reaching the threshold energy density eth hi 
the subsurface layer. P [t sp ik e > tfiow] is the probability 
of the spike lifetime being longer than the outflow time. 
Finally, P sp ut (E Q , Ri) is the probability of the cascade 
splitting into two cascades with each having not enough 
critical energy, or going too deep into the sample. Then, 

P(E ) =P(E ,R l ) = 

P [e (E ) > e th ] [1 - P S piit (E , R t )} P [t spike > t 

Note that at low values of E a , the threshold energy 
density should be e t h ~ 0.3eu o = 0.3nllo, and that the 
term P [t sp ik e > tfiow] will be roughly 1. This is because 
the crater is created mainly by "ballistic" events, not 
melt flow. On the other hand, at high energy density, the 
critical energy is much lower, e c ~ 3fcsT m ~ 0.07e[/ o , and 
spike times are important. From the MD simulations, 
we find that P sp ut (E a , Ri) depends on energy. It is near 
at low energies, and increases to 0.2 - 0.4 at higher 
energies. Of course, for energies much higher than here, 
it will eventually approach 1. 



a much simpler BCA calculation. The values for the 
calculated cratering probability compare well with the 
values from MD. 

For E a > 20 keV liquid flow is the main contribution 
and the energy density threshold for crater formation 
is lower, leading to an enhanced cratering probability. 
The cratering probability decreases rapidly for E > 50 
keV due to cascades being deep below the surface, and 
splitting due to fast recoils. Our simple model reproduces 
all these features, and extrapolating our results to E Q = 
200 keV gives P(200 keV) = 0.3 ± 0.1 and P(400 keV)= 
0.07 ± 0.05, which compares well with the experimental 
value of P(400 keV)= 0.03 from Ref. |1 



V. COMPARISON TO EXPERIMENT 

The crater radius from our MD simulations can be seen 
in Fig. p| as a function of bombarding energy. The radius 

1/3 

can be well approximated by a E dependence at low 
E Q , but at large E a there is a saturation, which is related 
to saturation of the energy deposition in the target. This 
energy deposition can be related to the known values 
for the stopping power, dE/dx, and it is included in the 
figure. There are some models that predict a dependence 
of the crater radius with (dE/dx) 1 ^ 2 (Refs. |50| , pT| ) and 
this is the dependence shown in Fig. || as R c = I Ax 1 / 2 I, 
where I = n -1 / 3 , and x = (dE / dxlIL/JJ) . 

The data from Birtcher et alJEiEl is also shown in 
Fig. ||. It is clear that within the uncertainties the 
experimental and simulated data agree very well, giving 
good confidence in the validity of our simulations. 

In Ref. |35| Donnelly and Birtcher discuss a criterion for 
crater formation based on the available energy density. 
They only look at high-energy cascades where spike times 
are much longer than the outflow time to form a large 
crater (with the mean size expected for those energies). 
This corresponds to our high energy region. 

The cratering efficiency is much smaller in the experi- 
ment (5%) than in our MD simulations (~ 50 %). There 
are two likely reasons for this: 



We show calculated values of P{E ) in Table III. 
The values in the table are calculated by using 
P [tspike > tfiow] = 1. P[e{E ) > e t h] and P sp u t were 
now obtained from MD, but could be obtained from 



1. Even though crystal orientation was not perfectly 
(4) known in the experiment, the incidence angle 6 was 

flow] about 15°, so there has been a significant amount 
of channeling (see Table || and note that both the 
mean range and straggling are larger for 8 = 15° 
than for 9 = 25° ), which decreases the energy 
deposition close to the surface and therefore the 
probability of cratering. 

2. The cratering probability of 0.05 was found from 
the ratio of cratering probability to crater destruc- 
tion cross section. This cross-section could be 
larger than estimated due to the large ion beam 
fluxes and possible enhanced atom mobility induced 
by the electron bombardment. 

Donnelly and Birtcher also assume thatr-the faceting 
of the crater sides occurs slowly by diffusion^. However, 
our simulations show (see Fig. 2) that the craters already 



are faceted directly in the collision cascades due to the 
crystal structure. 
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FIG. 4. (a) Total number of atoms above the surface, 
adatoms, and sputtered atoms as a function of the initial 
ion energy. Note that because many of the sputtered atoms 
sputter in hot clusters, and can subsequently evaporate from 
the cluster and redeposit on the surface, this initial sputtering 
yield does not exactly correspond to experimentally measured 
yields, (b) Crater rim length and rim width versus the initial 
ion energy. 



Expressing energy in eV and density in A 3 , A = 
(1.39 ±0.12) eV, B = (6.0 ±0.75) eVA 3 . Notice that 
B(2n/3)n = 0.74eV < A, The higher value of A 
indicates that the crater depth is typically smaller than 
the crater radius, as can be seen in Fig. || (a), where the 
open triangles represent the left side of Eq. H. 
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FIG. 5. (a) Number of adatoms, and total number of 
atoms outside the surface, vs. binding energy U. (b) Crater 
size vs. U. 



The number of total atoms above the surface after 
bombardment and the size of the crater rim can be seen 
in Figs. H (a) and (b). As expected from the behaviour of 
the crater radius, the number of atoms above the surface 
is linear with E Q at low energies and saturates at higher 
energies 

In the low energy regime we can cast our results in the 
following form: 



2tt\ n , E 



N cr = B 







E 
''TP 



(5) 



(6) 



VI. ROLE OF THE BINDING ENERGY OF THE 
TARGET 



Outside the ion beam community, it has generally 
been asspsaed, based on scaling laws, that crater size 
scales asElrE3 Uq 1 . This dependence was verified by 
macroscopic cratering events, like those of gas-gun 
experiments and_astronomical objects where gravity 
can be neglectecLj. However, this has been recently 
challenged by results in both EAM CuQ and LJ solidsE3, 
where the dependence was found to be Uq 2 . To test this, 
we repeated our 10 keV Xe bombardment simulations 
with the Au binding modified by changing the potential 



as explained in section II. 

In Fig. H a) the num ber o f atoms on top of the surface 
(as explained in section |ll D| ) is plotted versus the binding 
energy, and a quadratic dependence is found for bindings 
U in the range I.QUq-Q.SUq. This quadratic dependence 
may result from a combination of two factors. We found 
that the number of atoms in the melt and the lifetime of 
the spike scale as t/ _1 . Therefore, if the total number 
of adatoms scales as the product of the liquid atoms 
generated in the cascade and spike lifetime, the result 
depends on U as U~ 2 (See also Footnote |53|). In Fig. 
^-b the crater size is plotted as a function of the binding 
energy, also showing a dependence of U~ 2 . 



VII. COMPARISON WITH CLUSTER 
BOMBARDMENT 

Since there are a number of experiments and simula- 
tions dealing with cratering induced by cluster bombard- 
ment, we discuss some relevant cases. 

For cluster bombardment at low E c , the situation is 
such that the energy is deposited in a roughly hemi- 
spherical region and the threshold for crater formation 
is similar to what is found for atomic projectiles in the 
high energy regime. 

As mentioned in the introduction, Aderjan and 
UrbassekQ recently presented cratering results for 
Cu„ — >Cu where the number of atoms in the crater was 
found to scale as 



N cr = 131E(keV) - 656. 



(7) 



In addition, a scaling with U 2 was found. Using Eq. 
(0) for 20 keV bombardment of Cu, gives N cr = 1964, 
while using the density and binding of Cu in our Eq. 
(||) gives N cr = 1820 ± 230, in an excellent agreement 
with the previous estimate. Notice that the 20 keV Xe 
bombardment deposits all its energy close to the surface, 
as the Cu„ clusters do. The scaling found by Aderjan 
and Urbassek can not explain the experimental results of 
Xe bombardmenLof Au, nor the results for bombardment 
of Ceo on HOPGo, since the crater size increases always 
linearly with energy. This gives additional evidence for 
the need of a more complex model. Besides, Aderjan 
and Urbassek suggested that presence of viscosity may 
be the cause of the quadratic behavior with U . However, 
for EAM liquids near the freezing point, the viscosity 
v scales as v oc \JT m oc V~U (Ref. 55). If the crater is 
formed mainly by outflow we can consider some simplified 
cylindrical flow and use Poisson's equation, where the 
number N of flowing atoms is N tx V v~ x oc C/- 1 / 4 . This 
would give only an extra factor of 0.25 in the exponent 
and therefore it is quite unlikely to be the reason for the 
transition from linear to quadratic behavior. 

Typically, cluster bombardment simulations and ex- 
periments have been done at low energy per bombarding 
atom. For Cgo bombardment of HOPG, the crater 
"trace" (corresponding roughly to the rim length in our 
simulations), measured with an STM, was found to be 

1 /2 

proportional to [(dE/dx) n ] 1 , from 100 eV/atom up 
to 1 keV/atom, even though the experimental errors 



are quite larg E. A hemispherical crater is created, 
whose radius follows the law: (2/3) irR 3 nE c = E Q , with 
E c ~ 0.05£/o- This mean energy is consistent with the 
low energy densities found in the molten region. 

All these results support our findings that spikes play 
a major role in crater formation, and that a simple linear 
scaling of the crater volume with the bombarding energy 
is not a good description except at very low energies. 



VIII. COMPARISON WITH MODELS FOR 
OTHER KINDS OF MATERIALS 



Even though most cratering studies have been per- 
formed in the regime where the bombarding ion deposits 
most of its energy in_jelastic collisions, i-there are a 
number of experiments™] and simulationsEjE3 dealing 
with cratering in the regime where most of the energy 
is initially deposited in the electrons of the target. 
Experiments and simu.lst.ians of cratering in solids made 
of large biomoleculea£3'Ell suggested that the crater 
formation process is mainly due to the large pressure 
pulse following electronic relaxation. However, the. 
mechanisms of surface erosion on condensed gas solidsEa 
seem to be controlled by thermal spikes, with pressure 
playing only a secondary role in crater formation. 

Condensed gas solids and other soft materials can be 
reasonably well approximated by using simple 2-body 
potentials, as the Lennard- Jones potential (LJ). The 
track of excitations can be modeled as a cylindrical track 
of radius r cy u with atoms having some extra kinetic 
energy, i.e., a cylindrical spike. For tracks of fixed 
radius, r cy i ~ 21, with I = n -1 / 3 , there seems to be 
a critical energy density necessary for crater formation, 
which is close to nlMn. For rim formation, the energy 
density needed is even higher, of the order of the bulk 
modulus of the material or ~ 9nU. However, large 
craters can also be produced at a low energy density, 
as in the case of 100 keV Xe— > Au. In LJ solids this 
has been seen in simulations of tracks, where E c ~ 0.8/7 
produces no crater for a spike radius r sp ik e ~ 2Z, while 
for r S pike ~ 5Z the same energy density jean produce large 
craters because the spike lasts longerEa. LJ rare gases 
are comparatively stiffer than metals, and the pressure 
pulse associated with the high temperature spike takes a 
significant fraction of the energy when the energy density 
is large, but it can be neglected for low energy densities. 
This also decreases the relative lifetime of the spike as 
compared to metals. 

There are of course several differences to ion bom- 
bardment of metals, where the spike radius will vary 
significantly with the energy of projectile, and the 
energy deposition will not be uniform with depth, 
leading sometimes to energy deposition profiles closer to 
spherical geometry. In track simulations there is always 
a region at the surface, the top of the cylindrical track, 
which is energized. On the other hand, in the case of 
ion bombardment in the nuclear stopping regime, the 
peak in the energy deposition occurs below the surface. 
However, in the regime with energy above 10 keV for 
Xe bombardment, the created spike is initially roughly 



cylindrical, and one would hope certain scaling still to be 
valid. 

In EAM fee metals the critical energy density needed 
for crater formation for relatively narrow cascades is also 
close to nil , while for more extended cascades much lower 
energy densities also produce cratering. However, the 
energy density equivalent to the bulk modulus of Au is 
1.04 eV/A 3 = 4.54n£/ . This energy density is never 
reached in the Xe bombardment. On the other hand, the 
shear modulus is equivalent to 0.77nlV 

In the LJ track simulation, significant cascade splitting 
does not occur because of the initial conditions and can 
be neglected. In addition, because of the many body 
contribution in the EAM potential, the binding of small 
clusters is larger than in the LJ potential, and this may 
lead to an enhanced liquid flow to the surface. 

In Fig. [3] we compare the current results on cratering 
by Xe in EAM A^uwith those of cratering in a LJ solid by 
high-energy ionsEa. To make the two cases comparable, 
we give the abscissa in terms of (dE/dx)(Eo)(n~~ 1 / 3 /U 
for an incident ion with energy Eq, and the crater radius 
scaled by the characteristic length scale of the material 
n — 1/3_ rp ne gg ure snows that there is remarkably good 
agreement between the scaling of the crater radius R c 
with the stopping power. The rim width shown in the 
inset, however, does not follow any simple scaling law. 
The likely reason is the different heat spike geometry 
(cylindrical vs. hemispherical) leading to different crater 
shapes. 
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IX. CONCLUSIONS 



In conclusion, we have examined crater production 
by 0.4 - 100 keV Xe ions impacting on Au using 
classical MD simulations. On the basis of our results, 
we could show that the cratering mechanisms can not 
be understood in terms of a single parameter (such as 
a single energy density), but rather in terms of two 
energy regimes. In the low-energy regime (ion energies 
of ~ 1 — 20 keV in the present system) the cascades 
which produce craters need to have relatively high initial 
energy densities, ~ 0.5n[/o- However, at higher energies 
(> 100 keV) the liquid formed in the heat spike can 
become so long-lived that plenty of time is available for 
liquid flow, and a much lower initial energy density (~ 
0.2nUo) suffices for crater formation. These observations 
explains some apparent discrepancies between previously 
presented cascade models. 

We further demonstrated by direct simulation the 
importance of cascades splitting into subcascades on 
the cratering probability, and presented an analytical 
framework which can be used as a basis for further model 
development accounting for cascade splitting, energy 
density thresholds, and spike lifetimes. Furthermore, 
we showed that by using scaling laws and parameters 
calibrated from MD, simple BCA codes such as SRIM 
can be used to estimate cascade formation probabilities 

We obtained excellent agreement with crater sizes 
measured experimentally in the same system. 

Comparison with simulations of cluster bombardment 
and crater formation in LJ solids, which serve as models 
for some organic materials, showed that a wide range of 
microscopic cratering events can be understood in the 
same framework provided that appropriate scaling of the 
energy deposition and length scales are used. 

We also showed that macroscopic cratering laws be- 
have quite differently from the microscopic ones because 
of the importance of the liquid flow in the microscopic 
system, but not in the macroscopic systems. Conversely, 
we note that the results on cascade cratering may imply 
the need for a re-evaluation of macroscopic scaling laws 
for materials which expand strongly on melting, if the 
projectile can create a high enough energy density. 



X=(dE/dx)(n '7U) 

FIG. 6. Crater dimensions versus scaled dE/dx for both 
EAM Au and a LJ solid. The LJ results are from Ref. 



Thus we see that despite of the major differences 
between metals and some insulators which can be 
modeled as LJ solids, the crater sizes can still be 
understood in the same framework, when the incident 
energy deposited into atomic collisions (either directly or 
via electronic excitations) is considered along with the 
stopping power of the material. 
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